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A  COMPARISON  OF  EXPLICIT  TIME  INTEGRATION  TECHNIQUES 
FOR  THE  FINITE  ELEMENT  SHOCK  WAVE  EQUATIONS 


INTRODUCTION 

Many  physical  flow  processes  are  governed  by  nonlinear  hyperbolic  equations  of  the  form 

+  V  •  (VG)  -  F(V,  G)  U) 

6 t 

where  G  is  one  of  perhaps  several  conserved  quantities,  V  is  the  velocity  field  and  Fis  a  specified  func¬ 
tion  of  V  and  G.  Many  procedures  have  been  developed  for  obtaining  numerical  solutions  of  Eq.  (1). 
A  large  number  of  the  procedures  which  have  used  finite-difference  techniques  were  reviewed  by 
Roache  ( 1 )  in  his  classic  text  Computational  Fluid  Dynamics.  Other,  more  recent,  finite-difference  tech¬ 
niques  for  the  solution  have  been  reviewed  by  Sod  (2)  and  by  Book,  et.  al.  (3).  These  references 
clearly  show  that  very  accurate  and  very  sophisticated  finite-difference  techniques  have  been  developed 
for  the  solution  of  Eq.  (1).  However,  in  many  of  these  finite-difference  methods  a  serious  difficulty  is 
encountered  in  applying  the  methods  to  problems  with  complex  flow  boundaries.  This  can  be  particu¬ 
larly  true  for  the  more  accurate  methods  which  rely  on  the  use  of  a  "staggered"  mesh. 

A  more  natural  way  of  treating  complex  flow  boundaries  is  through  finite  element  discretization  of 
the  spatial  derivatives  in  Eq.  (1).  In  the  finite  element  method  —  a  very  good  introduction  is  given  in 
the  text  by  Baker  (4)  —  the  region  of  interest  is  divided  into  subspaces  (or  elements)  over  which  the 
dependent  variables  are  approximated  by  shape  functions.  The  constraints  on  the  shape  function 
coefficients  of  the  boundary  elements  then  naturally  provide  the  proper  treatment  of  complex  flow 
boundaries. 

Recent  work  by  many  researchers  —  see  for  example  The  Proceedings  of  the  First  International 
Conference  on  Numerical  Methods  in  Laminar  and  Turbulent  Flow  (5)  —  has  indicated  that  the  finite  ele¬ 
ment  method  (FEM)  has  considerable  promise  for  providing  numerical  solutions  for  fluid  flow  prob¬ 
lems.  However,  the  experience  base  in  using  the  FEM  in  fluid  flow  research  is  limited,  especially  when 
compared  with  finite-differences  and  there  are  as  yet  few  guideposis  to  help  a  researcher  select  the  most 
promising  paths.  The  need  for  systematic  evaluation  of  some  of  the  many  approaches  in  using  the 
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FEM  was  at  least  partially  met  by  the  work  of  Morrell,  et  al.  (6).  In  Ref.  (6),  the  solution  to  the  con¬ 
stant  velocity  advection  equation  was  obtained  by  the  use  of  four  finite  element  spatial  discretizations 
and  six  explicit,  two-step  time  integration  procedures.  Results  were  given  for  the  advection  of  two 
different  density  distributions;  the  first  a  cosine  hill  with  smooth  edges  and  the  second  a  square  hill  with 
abrupt  edges. 

In  the  present  work,  we  attempt  to  provide  further  guidance  by  a  systematic  evaluation  (much  like 
that  of  Ref.  6)  of  some  promising  approaches  for  solution  of  the  shock  wave  equations.  As  in  Ret  6, 
we  use  four  distinct  approaches  to  the  finite  element  spatial  discretization  and  six  procedures  for  the 
time  integration.  Since  the  shock  tube  problem  is  characterized  by  sharp  fronts  and  large  gradients,  the 
square  hill  results  of  Ref.  6  should  be  more  appropriate  for  comparison  with  the  present  results  than 
the  cosine  hill  results.  We  find  several  areas  of  substantial  agreement.  For  example,  the  Lax-W'endrofT 
second  order  accurate  time  integration  tends  to  yield  oscillatory  results.  Compared  with  the  Lax- 
Wendroff  time  integration,  the  Godunov  first  order  accurate  approach  is  quite  diffusive.  However,  here 
the  effect  of  the  diffusiveness  is  not  detrimental  but  is  beneficial  and  the  best  results  are  obtained  with 
the  Godunov  time  integration. 

NUMERICAL  METHOD 

The  governing  equations  of  gas  dynamics  are  given  in  many  different  forms  by  various  writers. 
For  example.  Bird,  Stewart  and  Lightfoot  (7)  give  them  as  below: 

continuity, 

'  ^-«-pV  V  Cal 

momentum  (neglecting  viscous  and  body  forces), 

nv 

pdT~~Vp  {2b) 

energy  (neglecting  viscous  forces  and  heat  flow  ). 


DE 

P~Di 


-  V  •  pX 


(2c) 


where  p  is  the  density,  V  is  the  velocity,  p  is  the  pressure  and  £  is  the  energy  per  unit  mass.  After 
some  rearrangement  of  terms,  Eqs.  (2)  may  be  rewritten  in  the  following  form  which  will  be  more 
appropriate  to  our  solution  procedure: 


+  V  •  pV  -  0 

Sr 

(3a) 

+  V  pVV-  -  Vp 
dt 

(3b) 

and 

+  V  •  eV  -  -  V  ■  p\ 
or 

(3c) 

where  eis  the  energy  per  unit  volume.  Equations  (3)  are  now  in  the  form  of  Eq.  (1). 

We  next  restrict  the  equations  to  one  dimension  and  shift  the  right  hand  side 

hand  side,  yielding 

term  to  the  left 

dt  dx 

(4a) 

+  f  W>  +  f  o»  -  o 

or  dx  dx 

(4b) 

+  —  (e«)  +  ~  (/>«)-  0 
dr  dx  dx 

(4c) 

where  u  is  the  velocity  in  the  x direction.  As  given  by  Sod  (2).  e  may  be  written  as 

f-pe  -  y  puJ  (5) 

where  «  is  the  internal  energy  per  unit  mass.  For  a  perfect  gas 

«  -  -Ei£-  (6) 

y  ~  l 

where  y  is  the  ratio  of  specific  heats.  Equations  (5)  and  (6)  may  be  combined  to  give  the  pressure 
explicitly  as 


P  “  (y  -  I)  e 


The  governing  equations  (4)  here  are  in  a  common  form  convenient  for  solution  as 

ia  +  i2.0 

dt  dx 


(7) 


(8) 


.  • 
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where 


P 

q  -  m  and  Q 
e 


m 

mu  +  p 
eu  +  pu 


and  where  m  —  pu.  We  now  turn  to  the  task  of  developing  an  appropriate  solution  procedure  based 
on  the  finite  element  method. 


Finite  Element  Solution  Pro*-  lure 

Much  as  was  done  in  Ref.  6,  we  shall  limit  the  scope  of  the  present  work  to  time  integration 
schemes  which  require  a  knowledge  of  the  dependent  variables  only  at  the  nodal  points  of  the  finite  ele¬ 
ment  grid.  Also,  we  will  consider  only  explicit  time  integration  schemes,  to  ensure  fast,  efficient  calcu¬ 
lations. 

In  solving  Eq.  (8)  by  the  finite  element  method,  the  domain  of  interest  is  divided  into  sub- 
domains  or  elements  of  finite  dimension.  The  dependent  variable  q  (here,  actually  variables)  are 
approximated  over  each  element  by  qe  as  follows 

q'ix.  t)  -  £rt;(r),V,(jc)  (9a) 

j- 1 

where  £f(r)  are  a  set  of  discrete  values  representing  the  dependent  variable  at  the  nodes  (or  grid 
points)  of  the  element,  and  A,(x)  (which  are  taken  to  be  of  the  same  form  for  each  element)  are  a  set 
of  interpolating  polynomials  (also  called  shape  or  basis  functions)  of  order  K  -  1  where  K  is  the 
number  of  nodes  per  element.  Likewise  the  variable  Q  (which  is  a  function  of  q  and  xand  t )  is  approx¬ 
imated  by  Qr  as 

Qe(x.  r)  -  £  T)(t).\,(x)  (9b> 

.i-i 

where  the  V.  are  the  same  interpolating  polynomials  as  in  Eq.  (9a)  and  the  7J  are  the  set  of  discrete 
values  of  Q'  at  the  nodal  points.  It  is  clear  that  since  qe  and  Qr  are  approximations  they  would  not 
satisfy  equation  (8)  but  would  give  some  error,  say  £r(x,  t).  Thus  for  the  approximations  on  each  ele¬ 
ment.  Eq  (8)  becomes 

~  (q‘)  +  ((?')  -  P(.v.  t).  (10) 

6/  dx 

Thus,  we  seek  a  solution  for  R;  and  7?  (for  any  assumed  q'  and  Q')  which  will  minimize  £'(.v.  /). 
The  minimization  of  F'ix,  i )  requires  that  the  error  £'(.v.  r )  be  orthogonal  with  A,(.v)  for  each  ele¬ 
ment,  therefore. 


JqL  E'(x,  t)Nj(x)dx  -  0  7  -  1 . K 


(11) 


where  L*  is  the  element  length.  Substitution  of  Eqs.  (9)  into  Eq.  (10)  and  then  into  Eq.  (11)  gives  the 


finite  element  discretization  equation  for  each  element,  i.e. 

(A/*‘)(/?'}  +  lKe]\Tr)  —  (0)  (12) 

where 

[AT]  -  m,j  -  J0L'  Nj(x)N,(x)dx  i,  j  -  I . K  (13a) 

[Ke]  -  ku  -  J0L'  Njlx).\',(x)dx  i.7-1 . K  (13b) 

{*')  -  f/?r,  . /?JC)r  03c) 

(r)  -(n.  T\ . 7Xlr  (13d) 


and  where  the  prime  denotes  spatial  differentiation  and  the  dot  denotes  time  differentiation.  When  Eq. 
(12)  is  written  for  each  element  and  when  the  element  connectivities  are  considered,  the  following  sys¬ 
tem  of  global  equations  are  obtained 

U/!(/?J  +  lAlID  -  10}  (14a) 

where 

R2 . A,v)r  (14b) 

{r|-{r,.  t2 . 7\lr  (i4c) 


and  where  l A/1  and  1  AT]  are  the  global  "mass"  and  "stiffness”  matrices  and  A  is  the  total  number  of 
nodes  in  the  discretization.  We  next  need  to  specify  the  forms  of  the  and  (A'l  matrices. 


5 


Linear  Elements 


One  of  the  simplest  elements  that  can  be  used  in  the  finite  element  method  is  the  linear  element; 
i.e..  a  first-order  interpolating  polynomial  and  for  Eqs.  (9)  K  ™  2  and  qe  and  Qr  will  be  given  by 

q'  -  JtfAf,  +  Re2S2 


and 


Qe  -  r,N,  +  r2.v2 


where 

.V,  -  1  -  (x/L"l 


-  (x/L'i  U  5) 

and  where  R\  and  Tf  are  the  values  of  qe  and  Qe  at  x  ■*  0  and  R2  and  7T  are  the  values  at  v  =  L‘  (Le 
is  the  element  length).  Substituting  Eq.  (15)  into  Eq.  (13a)  and  integrating  for  /,  j  “  1  to  2  gives  the 
elemental  mass  matrix 


(ATI 


U  2  1 
6  1  2 


(16) 


Likewise,  after  differentiating  Eqs.  (15)  and  substituting  into  Eq.  (13b).  the  elemental  "stiffness"  matrix 
is  obtained 


(*'] 


1  _1  1 

2  -1  1 


(17) 


Equation  (12)  can  then  be  written  for  each  element  as 


2  1 

//?! 

1 

-i  i 

77 

n 

1  2 

R'i 

+  T 

-l  i 

r: 

0 

The  assembled  mass  and  "stiffness"  matrices  for  the  linear  elements  are  given  in  Appendix  A 

Parabolic  Elements 

Also  commonly  used  is  the  parabolic  element,  i.e.,  one  using  a  second-order  interpolating  polyno¬ 
mial.  In  Eqs  (9)  K  —  3  and  qv  and  Qr  will  be  given  by 

q‘  -  Ky.V,  +  R'2\2+  RWy 
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and 


Q'  -  T\N\  +  T$N2  +  T\N 3 

where 

iV,  -  1  -  3(*/I')  +  2(x/Lr)2 

N2-  4(x/L')  -  4(x/Lf)2  (19) 

Ny-  -  ( x/Le )  +  2  (.<•///) 2 

and  where  R\,  R2,  R'i ,  Ff.  FJ,  and  T'S  are  the  va'ues  of  qe  and  Qe  at  x  —  0,  x  —  Le/2  and  ai  x  /. 

respectively.  Substitution  of  Eqs.  (19)  for  the  shape  functions  into  Eq.  (13a)  and  integrating  giv. 
following  for  the  elemental  mass  matrix 

4  2  -1] 

(A/0  -  2  16  2  .  (20) 

30  -1  2  4 

Likewise,  integration  of  Eq.  (13b)  after  differentiating  Eqs.  (19)  and  substitution  gives  the  following  for 
the  elemental  "stiffness"  matrix 

-3  4  -1 

«•!  -  -4  0  4  .  (211 

6  l  1  -4  ) 

Then  for  the  parabolic  element,  Eq.  (12)  can  be  written  as 

U_ 

30 

The  corresponding  assembled  mass  and  "stiffness"  matrices  are  given  in  Appendix  B. 

Condensed  IMI  Matrix  Formulation 

In  seeking  a  solution  of  Eqs.  (14a),  it  is  convenient  to  rewrite  it  as 

I*)--  lA/l-'U'llr}  (23) 

where  (A/l-1  denotes  the  inverse  of  (A/1.  Generally,  (A/)  is  a  large,  sparse  matrix  and  the  calculation 
of  (A/1'1  (which  is  a  full  matrix)  can  be  very  time  consuming  The  process  of  condensing  the  mass 


4  2-1 


1-3  4  -1 


2  16  2  R‘:  +  j  -4  0  4  F( 

1  2  4  1  -4  3  Ti 
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matrix,  i.e.,  collecting  the  off  diagonal  terms  to  the  diagonal,  is  a  popular  procedure  Ref.  (4)  for  avoid¬ 
ing  this  time  consuming  matrix  inversion. 


Thus,  we  define  the  condensed  mass  matrix  (A/,  ]  as  follows 


[Mc]  — 


v 

2>.,- 

i- 1 


mai  ~  0. 


/  -  1 . N 

i  ^  J\  i-  j  “  1 . /V 


(24) 


where  the  m„  are  the  entries  of  [A/1  and  A' is  the  rank  of  [A/].  The  matrix  [A/rl  now  replaces  1,4/1  in  Eq 
(14)  and  Eq.  (23)  becomes 

(/?}--  r).  (25) 

We  note  that  this  condensation  is  mass  preserving. 


With  either  Eq.  (23)  or  (25)  we  have  a  means  of  calculating  the  time  derivative  of  the  dependent 
variables  from  known  quantities.  What  remains  to  be  developed  are  the  procedures  for  advancing  the 
dependent  variables  in  time. 


Time  Integration  Schemes 

In  this  present  work  we  restrict  our  choice  of  methods  for  advancing  Ec  (23)  or  (25)  in  time  to 
standard  two  step  methods.  In  particular,  we  choose  two  methods  which  we  will  refer  to  as  the  Lav 
Wendroff  and  Godunov  methods  (see  Ref.  2). 

We  assume  that  the  nodal  variable  values  at  time  step  n  are  known,  i.e.,  (/?’’!  (and  thus  also 
(T"l).  The  solution  at  time  step  n  +  1  is  sought  and  the  time  increment  is  At.  The  standard  two  step 
method  can  be  written  as 

1st  step  I/r-°)  -  |K") +aA/(/r}  (26a) 

and 

2nd  step  iR"*'}  -  [R")  +  Ar{R"*"l  (26b> 

where  a  is  between  0  and  1.  The  time  derivative  {R"\  is  determined  from  Eq.  (23)  or  (25)  using  |R") 
and  thus  { 7~"J )  and  subsequently  [R'’+0|  is  determined  from  Eq  (23)  or  (25)  using  (R"*"|  (and  thus 
(r,*“))  calculated  from  Eq.  (26a).  In  the  solution  procedure  that  is  used,  we  also  considered  several 
types  of  weighting  or  smoothing  and  so  we  modify  Eq.  (26a)  by  premultiplying  |Ri  by  a  weighting 
matrix  [Ml,  and  Eq.  (26)  became 
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an 


1st  step  {/r+a|  -  [W]  {*"!  +  aAr[/T)  (27a) 

and 

2nd  step  {/Ti  -  !/r+,|  +  A?(/r*"|  (27b) 

and  we  see  that  if  (Hi  is  taken  to  be  the  identity  or  unit  matrix  (1),  then  Eqs.  (26)  and  (27)  are  exactly 
equivalent. 

The  two  well  known  methods  which  are  embodied  in  Eqs.  (27)  are  the  Lax-Wendroff  method 
when  a  “  1/2  and  the  Gudonov  method  when  a  —  1.  The  Lax-Wendroff  method  first  estimates  the 
time  derivative  at  the  mid-point  of  the  time  step  and  then  uses  the  mid-point  value  to  step  ahead  to  n 
+  1.  The  Goduov  method  first  estimates  the  time  derivative  at  a  full  time  step  ahead  and  then  uses 
this  estimate  to  advance  the  variabie  values  to  n  +  I. 

In  addition  to  the  standard  two-step  methods  represented  by  Eqs.  (26).  we  consider  the  two  types 
of  smoothing  in  the  time  integration  considered  in  Ref.  (6).  For  the  standard  weighting,  we  choose 
[1*1  “  Ul,  the  identity  matrix,  and  we  will  refer  to  the  standard  Lax-Wendroff  and  standard  Godunov 
methods. 

For  the  modified  two-step  methods  (to  use  the  terminology  of  Ref.  <6)).  the  weighting  matrix 
[  WJ  is  derived  from  the  global  mass  matrix  IA/J  and 

IM'J-'V-Tf  O-l . A’  (28a  l 


where 

C,  -  £  «„  C8b> 

j-  i 

and  where  the  m,;  are  the  terms  of  (A/1  and  A' is  the  rank  of  [.V/J.  Here  C,  is  the  sum  of  all  terms  in 
row  i  of  [A/1  and  thus  the  sum  across  row  i  of  l  is  unity.  For  the  modified  two-step  method.  ( H„.l 
replaces  [(*1  in  Eq.  (27a).  For  a  —  1/2  we  have  the  modified  Lax-Wendroff  method  and  for  o  «  1. 
the  modified  Godunov  method. 


The  second  type  of  weighting  that  we  consider  leads  to  the  "smoothed”  two-step  methods  For 
these  methods,  [Ml  -  [M'J  and 


in:'; 


M‘w<  “  0.  ;  —  1 . A' 

m„  i  j.  /.  ./  —  1 . X 

*  Ml  *  Q'  ' 
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where 


C,  -  £  mIJ  (29b) 

j-t 

j*' 

and  where  as  above  m,;  are  the  terms  of  [A/1  and  N  is  the  rank  of  [A/1  Here  C,  is  the  sum  of  the  terms 
of  the  i  th  row  of  [A/]  excepting  the  diagonal  term.  As  for  the  modified  and  standard  methods,  the  sum 
across  row  /'th  of  [W^l  is  unity.  In  the  smoothed  two-step  methods,  we  have  the  smoothed  Lax- 
WendrofT  method  -nd  the  smoothed  Godunov  method. 

Of  course,  for  each  of  the  two-step  methods,  we  consider  both  regular  and  condensed  finite- 
element  formulations  with  both  linear  and  parabolic  elements. 

TEST  PROBLEM 

As  a  test  problem  for  the  numerical  method,  we  have  selected  the  same  shock  tube  problem  and 
conditions  as  used  by  Sod  (2).  Initially  alt  velocities  (i.e.,  momenta)  are  zero,  the  pressure  and  density 
on  the  high  pressure  side  of  the  diaphragm  (region  1)  are  given  by  p\  “  p\  “  1.0  and  on  the  low  pres¬ 
sure  side  (region  5)  p5  “0.1  and  p5  —  0.125.  The  diaphragm  is  located  at  x  -  0.5  and  the  a- domain 
extends  from  0.0  to  1.0.  Except  for  some  special  test  cases  noted  below,  the  step  sizes  were  chosen  to 
be  Ax  -  0.01  and  A/  “  0.001,  yielding  a  Courant  number  of  0.22.  At  the  diaphragm  we  chose  to 
specify  the  initial  pressure  and  density  as  an  average  of  the  region  1  and  5  values,  i.e..  pd  0  “  0.55  and 
Pj  o  *  0  5625.  This  does  have  the  effect  of  initially  spreading  the  shock  over  two  zones  of  the  grid  but 
such  an  average  for  the  initial  conditions  seems  to  be  necessary  in  the  present  numerical  method. 
Numerous  tests  were  made  of  alternate  starting  conditions  including  variations  in  both  Ax  and  Ar.  The 
results  of  these  tests  will  be  discussed  later. 

At  each  end  of  the  shock  tube  we  chose  open  or  outflow  boundaries.  Since  our  concern  was  with 
the  propagation  of  the  shock  front  and  not  the  reflection  of  the  shock  from  a  closed  end  of  the  tube, 
these  boundaries  were  acceptable  for  this  work  and  slightly  simplify  the  problem.  A  small  difficulty  was 
encountered  in  some  of  the  calculations  in  the  form  of  small  numerical  oscillations  which  propagated  at 
about  two  times  the  shock  speed  and  could  reflect  off  the  open  end  of  the  tube  with  quite  adcerse 
effects.  In  order  to  ensure  that  a  smooth  solution  was  obtained  until  the  shock  passed  x  “  0.75.  4 
extra  points  were  added  at  each  end  of  the  tube  with  progressively  increasing  spacing  These  extra 
points  tended  to  damp  the  reflection  of  the  numerical  oscillations  from  the  open  ends  of  the  shock 
tube 
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NUMERICAL  RESULTS 


For  each  of  the  cases  considered,  we  show  the  present  numerical  results  plotted  against  the  ana¬ 
lytic  results  for  density,  velocity,  pressure  and  internal  energy.  The  analytic  solution  values  were  com¬ 
puted  for  the  grid  points  used  in  the  finite  element  formulation.  Thus,  even  for  the  analytic  solution 
the  shock  front  (and  also  the  contact  surface)  is  not  exactly  perpendicular  to  the  x-axis  but  is  plotted  as 
covering  one  zone.  We  categorize  the  present  results  according  to  the  type  of  weighting  used,  i.e., 
unweighted,  modified  weighting,  or  smoothed  weighting.  Within  each  category  we  discuss  the  results 
for  the  regular  finite  elements  method  (FEM)  and  the  condensed  finite  element  method  (CFM).  For 
each  of  these  we  consider  the  two  types  of  elements  and  the  two  types  of  time  integration. 

Unweighted  Method  Results 

For  the  results  in  this  section,  the  weighting  matrix  used  in  Eq.  (27a)  is  the  identity  or  unity 
matrix.  For  each  case  considered,  we  show  the  present  numerical  results  plotted  as  circles  and  the 
exact  analytic  solution  plotted  as  a  solid  line.  As  was  done  by  Sod  in  Ref.  (2)  we  show  the  density, 
velocity,  pressure  and  internal  energy  at  time  t“0.14  (i.e.,  at  the  time  when  the  shock  front  reaches  x 
-  0.75). 

The  first  results  that  we  consider  are  for  standard  FEM  for  the  Lax-Wendroff  (L  -  HI  type  of 
time  integration.  The  results  for  time  t=*=0.14  (i.e.,  when  the  shock  has  moved  to  x  —  0.75)  are 
shown  in  Fig.  1  for  the  linear  element  formulation  and  in  Fig.  2  for  the  parabolic  elements.  Several 
features  of  the  numerical  results  stand  out  in  these  figures.  The  most  obvious  are  the  large  oscilla¬ 
tions  in  the  solution  between  the  shock  front  and  the  foot  of  the  rarefaction  zone.  One  should  note 
that  (1)  this  method  (L  -  HT  is  non-diffusive,  (2)  there  is  no  weighting  or  smoothing  in  the  solution, 
and  (3)  there  is  no  artificial  diffusion  in  the  solution  procedure.  Here  the  shock  is  defined  over  only 
two  zones,  and  the  contact  surface  is  spread  only  over  two  zones.  Thus,  without  some  diffusivity.  these 
oscillations  are  not  surprising. 

Small  scale  oscillations  are  noted  at  about  x  —  0.1  and  0.9  These  are  the  numerical  oscillations 
mentioned  above  which  propagate  at  about  twice  the  shock  velocity  and  have  no  physical  significance. 
In  other  results  where  there  is  weighting  or  some  diffusivity  these  numerical  oscillations  are  either 
smaller  in  amplitude  or  absent  altogether.  We  also  note  that  the  numerical  results  agree  almost  exactly 
with  the  analytic  results  in  the  rarefaction  zone.  The  solution  with  the  parabolic  elements  shows  large 
oscillations,  a  larger  overshoot  in  front  of  the  shock  and  perhaps  also  a  steeper  shock  and  contact 
discontinuity  than  the  solution  with  the  linear  elements.  This  is  in  keeping  with  the  observation  of 
some  researchers  (e  g.  Baker,  Ref  8)  that  the  parabolic  element  formulation  seems  to  be  antidtlVusive 
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For  these  two  cases,  the  excellent  numerical  results  for  the  rarefaction  region  and  steepness  of  both  the 
shock  and  contact  surface  are  more  than  offset  by  the  large  oscillations  in  the  solutions. 

The  Godunov  time  integration  method  is  recognized  as  being  highly  diffusive  (see,  for  example, 
the  results  of  Morrell  in  Ref.  6).  The  numerical  results  for  the  standard  FEM  with  Godunov  time 
integration  are  shown  in  Figs.  3  and  4.  Compared  with  Figs.  1  and  2,  the  oscillations  are  nearly,  but 
not  quite,  eliminated.  We  see  several  additional  effects  of  the  diffusivity  of  the  method.  At  the  front 
of  the  rarefaction  zone,  at  x  —  0.35,  the  Godunov  results  do  not  as  accurately  follow  the  .talytic  solu¬ 
tion  as  do  the  L  -  W  results.  Also,  the  contact  discontinuity  is  spread  over  four  zones  instead  of  over 
two  zones.  The  results  for  the  linear  element  are  slightly  better  than  for  the  parabolic  element,  pri¬ 
marily  because  there  is  a  large  undershoot  in  the  internal  energy  at  x  ■*  0.76  for  the  parabolic  element, 
and  almost  no  undershoot  for  the  linear  element  case.  There  are  also  slightly  more  oscillations  in  the 
numerical  solution  for  the  parabolic  element.  Both  of  these  effects  appear  to  be  attributable  to  an 
antidiffusive  character  of  the  parabolic  element.  These  results  are  quite  good  and  we  will  compare 
these  results  with  some  of  the  results  in  Sod’s  paper  (Ref.  2)  after  we  discuss  other  results  of  the 
present  method. 
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Frg.  4  —  Finite  element  results  for  shock  tu be.  parabolic  elements.  Godunov  time  integration 


The  numerical  results  for  the  condensed  finite  element  method  (CFM)  are  shown  in  Figs.  5-8.  In 
comparison  with  Figs.  1-4  we  see  some  diffusivity  with  the  CFM  but  some  strong  adverse  effects  as 
well.  Numerical  results  are  not  shown  in  Fig.  6  because  of  the  unbounded  growth  of  the  spike  in  the 
solution  at  x  -  0.5.  Here,  the  antidiffusive  nature  of  the  parabolic  elements  combined  with  the  pecu¬ 
liarities  of  the  CFM  caused  the  solution  to  blow  up.  For  the  linear  elements  (Figs  5  and  7)  there  is  a 
small  spike  at  ,v  -  0.5  which  is  primarily  visible  in  the  plot  of  the  pressure.  For  the  parabolic  elements 
(see  Fig.  8)  the  spike  at  .v  —  0.5  is  quite  large  and  a  similar,  even  larger  spike  existed  in  the  numeric- 
solution  for  the  CFM  parabolic  element,  L  -  W  time  integration  at  t  =*  0.13  (ten  time  steps  earlier) 
The  Godunov  time  integration  did  slightly  damp  the  growth  of  the  spike  but  not  sufficiently,  as  that 
solution  blew  up  by  t  —  0.16  (twenty  steps  later).  For  the  linear  element,  the  shock  is  spread  over 
three  or  four  zones  and  the  contact  discontinuity  is  spread  over  six  or  seven  zones.  Considering  the 
apparent  diffusive  nature  of  the  CFM,  it  is  somewhat  surprising  that  the  linear  Godunov  case  has 
larger  overshoots  in  three  places  than  does  the  FEM  linear  Godunov  case,  (1)  behind  the  shock.  (2)  at 
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the  contact  discontinuity,  and  (3)  at  the  root  of  the  rarefaction  zone.  The  principal  benefit  that  might 
accrue  to  the  use  of  the  CFM  instead  of  the  standard  FEM  would  appear  to  be  if  an  adjustable  grid 
were  being  employed  and  the  [A/ 1  and  [A/I-1  matrices  needed  to  be  recomputed  frequently. 


Modified  Weighted  Method  Results 

The  two  principal  effects  of  the  use  of  the  modified  weighting  matrix  in  Eq.  (27a)  is  to  introduce 
some  damping  or  diffusi vity  into  the  system  and  to  more  strongly  connect  the  variable  values  at  each 
node  with  those  at  neighboring  nodes.  By  the  nature  of  the  mass  matrix,  this  later  effect  should  be 
more  pronounced  for  parabolic  elements  than  for  linear  elements.  Figures  9-12  show  the  numerical 
results  for  the  standard  FEM  when  the  modified  weighting  matrix  is  used. 

For  the  Lax-Wendroff  time  integration,  the  oscillations  have  been  greatly  reduced  but  are  still  not 
acceptable.  The  shock  is  spread  over  three  or  four  zones  and  the  contact  discontinuity  is  spread  over 
five  or  six  zones.  A  mild  overshoot  exists  for  both  elements  at  x  —  0.5.  These  results  seem  to  suggest 
that  some  artificial  viscosity  for  the  L  —  ^integration  might  be  very  useful. 

For  the  Godunov  time  integration,  the  coupling  between  neighboring  points  leads  to  some 
overshoots  (at  x  =  0.5  and  0.7)  although  for  the  parabolic  element,  these  are  not  at  all  severe.  The 
shock  is  spread  over  two  or  three  zones  and  the  contact  discontinuity  is  spread  over  four  or  six  zones. 
In  all  four  cases,  the  damping  due  to  the  modified  weighting  prevents  the  numerical  results  from  fol¬ 
lowing  the  corner  at  the  rarefaction  front  at  x  *  0.3  as  well  as  was  done  for  the  unweighted,  standard 
FEM.  As  was  true  in  the  results  discussed  above,  the  shock  front  and  contact  discontinuity  were  spread 
over  fewer  elements  when  the  parabolic  element  was  used  than  where  the  linear  element  was  used.  In 
particular,  the  plateau  in  the  internal  energy  between  x  —  0.62  and  0.74  is  significantly  better  defined 
with  the  parabolic  element. 

For  the  modified  weighting  matrix  used  with  the  condensed  finite  element  formulation,  damping 
is  again  evident  in  the  results,  which  are  shown  in  Figs.  13-16.  The  numerical  results  seem  to  be  con¬ 
sistently  poorer  than  for  the  standard  FEM.  For  the  linear  elements  (Figs.  13  and  15)  the  contact 
discontinuity  is  nearly  obscured  (see  the  plots  of  density)  and  except  for  the  parabolic  element. 
Godunov  time  integration  case,  (he  plateau  in  the  internal  energy  behind  the  shock  is  hardly  resolved  at 
all. 
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Smoothed  Weighting  Method  Results 

The  numerical  results  for  the  use  of  the  smoothed  weighting  matrix  in  Eq.  (27a)  are  shown  in 
Figs.  17-20.  Morrell  (6)  had  noted  that  at  times  this  method  seemed  to  be  unstable.  In  the  present 
results  a  large  undershoot  occurs  in  the  density  and  pressure  at  x  ■>  0.5  and  is  sufficiently  severe  to 
blow  up  the  four  runs  with  parabolic  elements  considerably  before  the  time  when  the  shock  reached  x 
**  0.75.  Thus,  only  results  for  linear  elements  are  shown  here.  Even  though  this  technique  seems  to 
be  the  most  diffusive  of  all  of  the  methods  considered  (note  the  rarefaction  front  at  x  «  0.3),  the 
results  have  some  of  the  larger  overshoots  (or  undershoots).  The  undershoot  at  x  **  0.5  is  particularly 
bad  and  was  noted  above.  The  overshoots  behind  the  shock  are  also  quite  bad  (note  the  velocity  plots). 
In  each  of  the  velocity  plots  one  or  two  points  (at  x  =  0.51  and  0.52)  are  off  the  graph.  These  points 
coincide  with  the  undershoots  in  the  density,  pressure  and  internal  energy.  Considering  the  present 
results  along  with  the  results  in  Ref.  6,  this  weighting  technique  seems  to  have  few  advantages  to 
recommend  it. 

Evaluation  of  Results 

The  evaluation  of  results  relies  primarily  on  a  visual  comparison  of  the  numerical  and  analytic 
results.  Primarily  because  of  the  strong  oscillations,  none  of  the  Lax-Wendroff  time  integration  cases 
are  included  in  our  group  of  belter  results.  Thus  all  those  discussed  below  are  with  Godunov  time 
integration.  The  best  two  cases  are  for  standard  FEM  and  standard  no  weighting  or  averaging  as 
shown  in  figures  3  and  4.  By  a  small  margin  the  linear  element  (Fig.  3)  gives  the  better  results.  The 
next  two  best  cases  are  for  standard  FEM  with  the  modified  weighting  matrix  as  in  Figs.  11  and  12. 
Mere,  the  parabolic  element  gives  slightly  better  results  than  the  linear  element.  The  next  best  results 
are  for  the  CFM,  modified  weighting  matrix  with  parabolic  element.  Three  other  cases  are  sufficiently 
good  to  merit  a  fair  rating.  These  are  the  ones  using  the  CFM  formulation  and  linear  elements  and  <  1 ) 
no  averaging  (Fig  7),  (2)  modified  weighting  (Fig.  15)  and  (3)  the  smoothed  averaging  (Fig.  20). 

For  this  class  of  problem,  contrasted  with  the  simple  advection  problem  considered  by  Morrell 
(Ref.  6),  the  Godunov  time  integration  is  clearly  superior  to  the  Lax-Wendroff  time  integration.  The 
use  of  artificial  viscosity  might  well  improve  the  results  for  the  L  —  W  approach  but  would  most  likely 
be  detrimental  at  the  rarefaction  front.  The  challenge  would  be  to  obtain  results  as  good  as  those  in 
Figs.  3,  4  or  12  Overall,  the  linear  and  parabolic  elements  seem  to  perform  almost  equally  well,  since 
often  the  antidiffusive  nature  of  the  parabolic  element  needs  to  be  offset  with  some  artificial  diffusion 
or  damping  The  results  of  the  condensed  mass  matrix  formulation  suggest  that  this  approach  might  be 
useful  under  some  circumstances,  such  as  when  an  adaptive  grid  is  used  and  a  new  matrix  inverse 
would  need  to  be  computed  each  time  the  finite  element  grid  is  altered.  However,  for  calculation  with 
a  fixed  grid  the  CFM  approach  is  not  recommended. 

22 


PRESSURE  DENSITY  PRESSURE  DENSITY 


Comparison  with  Finite-Difference  Results 

While  we  do  not  show  a  direct  comparison  with  results  of  finite-difference  methods,  a  quantitative 
comparison  is  still  possible  and  reasonable.  For  this  comparison  we  shall  consider  the  results  of  the 
best  four  cases  of  the  present  methods  and  the  results  of  the  best  four  of  the  cases  presented  in  Ref  2 
by  Sod.  Sod  notes  that  the  best  two  without  corrective  procedures  are  the  ones  for  Godunov's  and 
Hyman’s  methods.  Hyman’s  method  spreads  the  shock  over  four  zones  and  the  contact  discontinuity 
over  seven  or  eight.  Godunov’s  method  requires  five  or  six  zones  and  the  contact  discontinuity  seven 
or  eight.  Consequently,  for  both  finite-difference  methods,  the  flat  crest  in  the  internal  energy  behind 
the  shock  is  poorly  defined  and  quite  rounded.  We  would  compare  these  FD  results  with  our  results  in 
Figs.  11  and  12.  Our  results  show  the  shock  over  three  zones  and  the  contact  discontinuity  over  four 
zones  for  the  parabolic  element  (two  elements)  and  over  six  or  seven  elements  for  the  linear  case. 
While  our  results  do  have  some  small  oscillations  and  some  slight  overshoots,  it  would  seem  that  these 
are  a  fair  trade  for  the  crispness  in  the  definition  of  the  shock  and  contact  discontinuity. 

Of  the  corrective  procedures  which  he  considered,  it  seems  that  Sod  prefers  the  artificial  compres¬ 
sion  method  (ACM)  to  the  antidiffusion  method.  It  appears,  however,  that  in  the  results  he  presents, 
that  the  antidiffusion  method  (used  with  the  two  step  Lax-Wendroff)  gives  better  results  than  either  the 
hybrid  method  with  ACM  or  Godunov’s  method  with  ACM.  The  antidiffusion  method  defines  the 
shock  over  two  zones  and  the  contact  discontinuity  over  six  or  seven  zones.  Also  there  exists  a  slight 
overshoot  at  the  right  corner  of  the  rarefaction  zone.  We  would  compare  these  results  with  our  results 
for  Godunov  time  integration  using  standard,  unweighted  finite  elements  —  particularly  with  the  linear 
element.  There  are  slight  oscillations  in  the  FEM  results  (there  being  no  artificial  diffusion  in  the 
present  numerical  results  as  there  was  in  the  Ref.  2  numerical  results)  and  a  very  small  overshoot  in 
front  of  the  shock.  However,  the  shock  is  resolved  over  two  zones  and  the  contact  discontinuity  is 
spread  only  over  four  zones.  Thus  even  with  the  oscillations  in  the  internal  energy  behind  the  shock, 
the  flat  peak  is  somewhat  better  defined  than  it  was  in  the  Ref.  2  results  for  the  antidiffusion  method. 

We  certainly  recognize  that  there  have  been  many  improvements  in  finite-difference  methods 
since  Sod's  1978  review  article,  and,  with  these  improved  finite-difference  methods,  better  results  can 
be  obtained  for  this  shock  tube  problem.  The  present  results  obtained  using  the  finite  element  method 
seems  to  be  very  good,  especially  considering  that  no  corrective  procedures  were  employed. 

Tests  of  Alternate  Initial  Conditions 

In  the  earlier  discussion  of  the  initial  conditions,  we  noted  that  the  initial  pressure  and  density 
values  were  averaged  at  the  diaphragm  location.  For  example,  at  .v  ~  0.49,  0.50  and  0.51,  the  values  of 


25 


pressure  were  p  “  1.0,  0.55  and  0.1.  The  density  was  similarly  averaged.  We  have  seen,  that  for 
these  initial  conditions,  the  computer  program  runs  quite  well.  We  also  tried  initial  conditions  without 
averaging  at  the  diaphragm  location,  i.e.,  at  x  —  0.49  and  0.50,  p  —  1.0  and  0.1  and  p  *  1.0  and  0.125. 
These  cases  did  not  run  well  at  all.  This  is  somewhat  puzzling,  since  the  latter  seem  to  be  the  initial 
conditions  most  other  researchers  employ.  In  order  to  understand  this  situation  better  we  ran  a  series 
of  additional  cases. 

first  we  ran  four  cases  for  the  standard  finite  element  formulation  with  linear  elements  and 
Godunov  time  integrations  and  with  pressure  and  density  averaging  at  x  —  0.50.  With  Ax  —  0.02 
(instead  of  Ax  —  0.01),  the  program  ran  quite  well,  but  the  oscillations  in  the  solution  were  slightly 
larger  and  the  shock  and  contact  discontinuity  were  spread  over  a  slightly  greater  ,v  distance.  With  Ax 
—  0.005  (uniformly),  the  program  again  ran  quite  well  and  the  solution  was  slightly  improved  over  the 
standard  case  of  Ax  —  0.01.  The  oscillations  were  smaller  and  the  shock  and  contact  discontinuity  were 
portrayed  more  steeply. 

We  also  ran  two  cases  with  the  standard  Ax  ~  0.01  but  with  one  or  four  extra  points  added  at  or 
near  the  diaphragm;  at  x  ■  0.505  (and  0.485,  0.495  and  0.515).  In  each  of  these  two  cases  the 
diaphragm  was  moved  from  x  —  0.50  to  x  —  0.505  and  the  density  and  pressure  values  were  averaged 
at  x  »  0.505.  This  gave  a  steeper  initial  gradient  than  standard  but  with  the  midpoint  still  defined.  The 
two  cases  ran  quite  well  matching  the  results  in  Fig.  3  except  for  having  smaller  oscillations  a  x  “  0.51. 
For  these  four  cases,  the  initial  gradient  in  pand  p  differed  by  a  factor  of  four  with  little  difference  in 
the  numerical  results  other  than  the  solution  being  smoother  with  a  closer  grid  spacing 

In  the  tests  without  averaging  of  p  and  p  at  the  diaphragm  location,  we  used  the  standard  finite 
element  formulation,  the  condensed  mass  matrix  formulation,  no  smoothing,  modified  smoothing. 
Lax-Wendro ff  and  Godunov  time  integration,  linear  and  parabolic  elements,  and  different  values  of  A.v 
and  At.  None  of  these  test  cases  would  provide  a  solution  beyond  t  -  0.05.  In  each  case,  a  severe 
undershoot  occured  in  the  density  and  pressure  at  the  diaphragm  location  (x  “  0.5)  as  the  solution 
developed.  On  each  side  of  x  —  0  5  the  values  of  pand  p  would  begin  to  change  to  meet  those  from 
the  other  side  but  at  x  —  0.5  the  values  would  plunge  to  zero.  One  of  the  longer  running  cases  was 
with  Ax  —  0.02  (giving  an  initial  gradient  matching  the  standard,  averaged  case  with  A.v  —  0.01).  This 
case  was  for  standard  FEM,  linear  element,  Godunov  time  integration  and  no  smoothing  and  ran 
beyond  t  “  0.04  but  not  to  t  “  0.05.  Other  cases  that  ran  as  well  were  in  a  group  with  A.v  —  0.01  (the 
standard  Ax  values),  standard  FEM  and  with  the  modified  weighting  matrix.  This  latter  group  was  run 
with  both  Godunov  and  Lax-Wendroff  time  integration  techniques  and  with  A t  ”  0.001,  0.005.  and 
0.0001.  The  smaller  values  of  A/ seemed  to  help  but  only  by  a  small  amount  The  use  of  the  damping 


in  the  modified  weighting  matrix  clearly  helped,  but  it  did  not  help  enough  for  the  computer  program 
to  run  to  completion. 

ADDITIONAL  NUMERICAL  RESULTS 


The  results  presented  in  the  previous  sections  were  all  for  low  pressure  and  density  ratios  across 
the  diaphragm  (10:1  or  less).  On  the  high  pressure  side,  pand  p  were  given  by  ps  and  p5  ■»  1.0,  while 
on  the  (ow  pressure  side  of  the  diaphragm  the  initial  conditions  were  p\  »  0.1  and  pi  ~  0.125.  These 
conditions  were  chosen  to  match  the  initial  conditions  used  by  Sod  (2)  in  his  paper,  and  were  appropri¬ 
ate  conditions  for  evaluating  the  various  integration  methods.  However,  these  conditions  did  not 
severely  test  the  numerical  method.  In  order  to  evaluate  the  method  more  fully,  additional  calculations 
were  made  with  higher  density  and  pressure  ratios.  These  additional  calculations  were  made  with  linear 
elements,  Godunov  time  step  integration  and  no  weighting. 


For  these  additional  numerical  results,  calculations  were  made  with  pressure  ratios  as  high  as  10’:1 
and  with  density  ratios  as  high  as  500:1.  Results  from  some  of  these  calculations  are  shown  in  Figs. 
21-30.  There  were  no  significant  problems  encountered  in  making  the  calculations  for  the  higher  pres¬ 
sure  ratios,  but  problems  were  encountered  in  trying  to  make  calculations  with  higher  density  ratios.  It 
appeared  that  some  damping,  in  addition  to  that  provided  by  the  Godunov  integration  method,  would 
be  needed  in  order  to  obtain  successful  calculations  for  higher  density  ratios. 


x  -  distance  x  -  distance 

Eig  21  —  lliph  pressure  resulls  for  shock  luhe.  linear  elements,  (ioilunm  lime  inlegunon.  no  wcighling. 
p,  -  Pl  -  I  0.  Ps  -  10-’.  p,  -  100 1) 
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Fig  30  —  High  pressure  results  for  shock  tube,  linear  elements.  Godunov  time  integration,  no  weighting. 

Pi  -  Pi  -  300  0.  ps  -  105.  ps  -  3000 

Jn  the  following  figures,  the  density  and  pressure  on  the  low  pressure  side  of  the  diaphragm  were 
set  to  1.0.  On  the  high  pressure  side,  the  densities  were  either  10  or  100  and  the  pressures  were  varied 
from  100  to  10'  Figures  21-26  show  the  results  for  p5  -  10:  to  107  with  p<  —  100  There  are  several 
principal  features  in  these  figures.  First,  there  is  an  increased  oscillation  in  front  of  the  shock  front  as 
the  pressure  is  increased.  There  is  also  a  decreased  distance  between  the  shock  front  and  the  fool  of 
the  rarefaction  zone.  Consequently,  the  features  in  this  region  are  not  as  well  defined  as  for  the  low 
pressure  case. 

Somewhat  similar  results  are  given  in  Figs.  27-29  for  p<  —  10  and  p<  ■*  10:,  103  and  10'.  In  all  of 
these  cases,  the  shock  is  defined  quite  well  by  the  results  of  the  numerical  method,  but  with  some  lead¬ 
ing  oscillation,  and  the  numerical  and  analytical  results  agree  quite  well  in  the  rarefaction  zone  The 
contact  discontinuity  is  spread  over  four  of  the  linear  elements  which  is  comparable  to  the  low  pressure 
results 

Figure  30  shows  the  numerical  results  for  a  case  near  the  limits  of  what  the  numerical  method 
seems  to  be  capable  without  the  addition  near  the  shock  of  artificial  viscosity  or  other  additional  damp¬ 
ing  In  this  case,  p<.  and  p>  -  300  Calculations  were  successfully  made  for  some  higher  density  ratios 


32 


but  only  with  a  corresponding  reduction  in  the  pressure  ratio.  The  oscillations  near  the  shock  are 
stronger  than  for  the  lower  density  cases  as  are  the  oscillations  at  the  foot  of  the  rarefaction  zone.  At 
higher  densities,  the  oscillations  at  the  shock  rapidly  grow  and  overwhelm  the  entire  solution.  For 
these  present  conditions,  however,  the  solution  is  still  reasonably  well  behaved,  and  the  numerical  and 
analytical  results  agree  rather  well. 

The  purpose  of  these  additional  cases  was  not  to  define  the  full  range  of  initial  conditions  to 
which  the  numerical  method  is  applicable  but  rather  to  determine  if  the  method  could  be  applied  to  a 
broader  range  of  conditions  than  those  considered  by  Sod  (2).  The  results  shown  in  Figs.  21-30  do 
show  that  the  numerical  method  is  capable  of  application  to  a  wide  range  of  initial  conditions  and  that  it 
does  give  results  which  are  in  good  agreement  with  analytical  results. 

CONCLUSIONS 

While  the  present  work  certainly  does  not  indicate  a  best  approach  in  applying  the  finite  element 
method,  valuable  insights  have  been  achieved.  First,  the  FEM  is  capable  of  providing  very  good  results 
for  fluid  flow  problems  such  as  the  shock  tube.  While  parabolic  elements  have  a  potentially  useful 
antidiffusive  characteristic,  they  must  be  used  with  care  since  often  the  antidiffusive  nature  needs  to  be 
balanced  by  some  damping.  The  linear  element  is  to  be  recommended  for  its  simplicity  and  its  lack  of 
either  diffusivity  or  antidiffusiveness.  The  modified  weighting  or  averaging  approach  has  led  to  some 
very  good  results  here,  but  one  might  well  prefer  the  addition  of  specific  artificial  diffusion,  even  if  the 
amount  of  artificial  diffusion  is  problem  dependent.  The  condensed  mass  matrix  formulation  of  the 
finite  element  method  seems  to  have  both  a  diffusive  as  well  as  a  compressive  nature.  It  might  be 
recommended  if  many  inversions  of  the  "mass”  matrix  should  be  needed,  but  otherwise  the  CFM 
approach  should  be  used  with  considerable  caution. 
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Appendix  A 

LINEAR  ELEMENT  MATRICES 


Let  the  x  domain  be  subdivided  into  N  linear  finite  elements  each  of  tength  U.  The  nodes  are 
numbered  consecutively  from  1  to  N  +  1  as  x  increases  from  0  to  L.  The  element  mass  [Mf  1  and 
advection  [Ke  1  matrices  are  given  by  Equations  (16)  and  (17). 

The  assembled  mass  [A/]  and  advection  [A'l  matrices  of  the  N  degree  of  freedom  system  are 


found  as 
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Appendix  B 

PARABOLIC  ELEMENT  MATRICES 


Let  the  x  domain  be  subdivided  into  N  parabolic  finite  elements  each  of  length  L'.  The  nodes  are 
numbered  consecutively  from  I  to  2A'  +  1  as  x  increases  from  0  to  L.  The  even  numbered  nodes 
correspond  to  element  mid-point  nodes.  The  element  mass  lA/‘l  and  advection  [A'')  matrices  are  given 
by  Equations  (20)  and  (21). 

The  assembed  mass  [A/]  and  advection  (Al  matrices  of  the  2, V  degree  of  freedom  system  are 


found  as 


